Advanced SPPA Methods - NCKDE

Author

Federico Jose Rodriguez

Published

September 9, 2024

Modified

September 9, 2024

Getting Started

Data Sources

Data for this exercise are from public sources and will be used to analyse the distribution of childcare centers in the Punggol planning area. Two datasets in ESRI shapefile format will be used:

  • A line feature geospatial dataset which includes the road network of Punggol planning area

  • A point feature geospatial dataset which includes the location of childcare centers in the Punggol planning area

Installing and launching R packages

This exercise will make use of five R packages: sf, spNetwork, tidyverse, and tmap.

  • sf - for importing, managing and processing vector-based geospatial data

  • tidyverse - collection of packages for performing data importation, wrangling and visualization

  • tmap - for plotting cartographic quality maps

  • sPNetwork - provides functions for performing SPPA methods like KDE and K-function on a network. The package can also be used to build spatial matrices to conduct traditional spatial analyses with spatial weights based on reticular distances

The code chunk below uses p_load() of pacman package to check if the packages are installed in the computer. It installs them first if they are not. It then loads them into R.

pacman::p_load(sf, spNetwork, tmap, tidyverse)

We use the random seed 1234 to ensure reproducibility of results

set.seed(1234)

Data Import and Preparation

The code chunks below uses st_read() of the sf package to load the street and childcare data into their respective dataframes.

network <- st_read(dsn="data/geospatial", 
                   layer="Punggol_St")
Reading layer `Punggol_St' from data source 
  `C:\drkrodriguez\ISSS626-GAA\In-class\In-class_Ex03\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 2642 features and 2 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: 34038.56 ymin: 40941.11 xmax: 38882.85 ymax: 44801.27
Projected CRS: SVY21 / Singapore TM
childcare <- st_read(dsn="data/geospatial",
                     layer="Punggol_CC")
Reading layer `Punggol_CC' from data source 
  `C:\drkrodriguez\ISSS626-GAA\In-class\In-class_Ex03\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 61 features and 1 field
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 34423.98 ymin: 41503.6 xmax: 37619.47 ymax: 44685.77
z_range:       zmin: 0 zmax: 0
Projected CRS: SVY21 / Singapore TM

We examine the contents of these objects by running the code chunks below

network
Simple feature collection with 2642 features and 2 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: 34038.56 ymin: 40941.11 xmax: 38882.85 ymax: 44801.27
Projected CRS: SVY21 / Singapore TM
First 10 features:
     LINK_ID                   ST_NAME                       geometry
1  116130894                PUNGGOL RD LINESTRING (36546.89 44574....
2  116130897 PONGGOL TWENTY-FOURTH AVE LINESTRING (36546.89 44574....
3  116130901   PONGGOL SEVENTEENTH AVE LINESTRING (36012.73 44154....
4  116130902   PONGGOL SEVENTEENTH AVE LINESTRING (36062.81 44197....
5  116130907           PUNGGOL CENTRAL LINESTRING (36131.85 42755....
6  116130908                PUNGGOL RD LINESTRING (36112.93 42752....
7  116130909           PUNGGOL CENTRAL LINESTRING (36127.4 42744.5...
8  116130910               PUNGGOL FLD LINESTRING (35994.98 42428....
9  116130911               PUNGGOL FLD LINESTRING (35984.97 42407....
10 116130912            EDGEFIELD PLNS LINESTRING (36200.87 42219....
childcare
Simple feature collection with 61 features and 1 field
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 34423.98 ymin: 41503.6 xmax: 37619.47 ymax: 44685.77
z_range:       zmin: 0 zmax: 0
Projected CRS: SVY21 / Singapore TM
First 10 features:
      Name                      geometry
1   kml_10 POINT Z (36173.81 42550.33 0)
2   kml_99 POINT Z (36479.56 42405.21 0)
3  kml_100 POINT Z (36618.72 41989.13 0)
4  kml_101 POINT Z (36285.37 42261.42 0)
5  kml_122  POINT Z (35414.54 42625.1 0)
6  kml_161 POINT Z (36545.16 42580.09 0)
7  kml_172 POINT Z (35289.44 44083.57 0)
8  kml_188 POINT Z (36520.56 42844.74 0)
9  kml_205  POINT Z (36924.01 41503.6 0)
10 kml_222 POINT Z (37141.76 42326.36 0)

We see from the output of calling childcare that the dataset includes a third dimension– a z-coordinate. The spNetwork package can only accept objects in (x,y) and not in (x,y,z)

To go around this, we re-import the data with st_zm() (using default arguments: drop = TRUE, what = “ZM”) to remove the z-coordinate

childcare <- st_read(dsn="data/geospatial",
                     layer="Punggol_CC") %>%
  st_zm()
Reading layer `Punggol_CC' from data source 
  `C:\drkrodriguez\ISSS626-GAA\In-class\In-class_Ex03\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 61 features and 1 field
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 34423.98 ymin: 41503.6 xmax: 37619.47 ymax: 44685.77
z_range:       zmin: 0 zmax: 0
Projected CRS: SVY21 / Singapore TM

Checking the content of the object shows that the third dimension is now dropped

childcare
Simple feature collection with 61 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 34423.98 ymin: 41503.6 xmax: 37619.47 ymax: 44685.77
Projected CRS: SVY21 / Singapore TM
First 10 features:
      Name                  geometry
1   kml_10 POINT (36173.81 42550.33)
2   kml_99 POINT (36479.56 42405.21)
3  kml_100 POINT (36618.72 41989.13)
4  kml_101 POINT (36285.37 42261.42)
5  kml_122  POINT (35414.54 42625.1)
6  kml_161 POINT (36545.16 42580.09)
7  kml_172 POINT (35289.44 44083.57)
8  kml_188 POINT (36520.56 42844.74)
9  kml_205  POINT (36924.01 41503.6)
10 kml_222 POINT (37141.76 42326.36)

Visualization of the sf objects

The code below uses plot() in one map. The add=T argument in the second line allows the two plots to be added one over the other.

plot(st_geometry(network))
plot(childcare, add=T, col="red", pch=10)

The following code chunk produces a similar map using tmap package.

tmap_mode('view')
tmap mode set to interactive viewing
tm_shape(network) +
  tm_lines() +
  tm_shape(childcare) +
  tm_dots(col = "red")
tmap_mode('plot')
tmap mode set to plotting

Network Constrained Spatial Point Analysis

Preparing the lixels objects

A requirement for NKDE is that the lines objects needs to be cut into lixels. The code below uses lixelize_lines() on network, using a length of 700 and a minimum distance (mindist) of 350.

lixels <- lixelize_lines(network,
                         700,
                         mindist = 350)

Generating the line center points

The code below generates the center points for the lixels

samples <- lines_center(lixels)

xxx

tmap_mode("view")
tmap mode set to interactive viewing
tm_shape(lixels) +
  tm_lines() +
  tm_shape(samples) +
  tm_dots(size = 0.01)
tmap_mode("plot")
tmap mode set to plotting

Computing NKDE

The code below computes for the NKDE of childcare centres around network

densities <- nkde(network, 
                  events = childcare,
                  w = rep(1, nrow(childcare)),
                  samples = samples,
                  kernel_name = "quartic",
                  bw = 300, 
                  div= "bw", 
                  method = "simple", 
                  digits = 1, 
                  tol = 1,
                  grid_shape = c(1,1), 
                  max_depth = 8,
                  agg = 5, 
                  sparse = TRUE,
                  verbose = FALSE)

The output of the code is a list of density or intensity values. We copy these values into the original lixel and lixel midpoint dataframes.

samples$density <- densities
lixels$density <- densities

Note that the previous values are based on metres and resulted in very low density values. We can change the density values to event per square km by using the code below

samples$density <- samples$density*1000
lixels$density <- lixels$density*1000

Computing the K and G Functions

The code block below computes for the K- and G-functions based on the data using kfunctions()

kfun_childcare <- kfunctions(network, 
                             childcare,
                             start = 0, 
                             end = 1000, 
                             step = 50, 
                             width = 50, 
                             nsim = 49, 
                             resolution = 50,
                             verbose = FALSE, 
                             conf_int = 0.05)

The output can be visualized by calling plotk for the K-function and plotg for the G-function

kfun_childcare$plotk